v20200329.2025
Os exemplos aqui apresentados estão disponíveis. Caso queira usá-los, crie um projeto, coloque o arquivo desta aula e os seguintes arquivos na pasta do mesmo:
É uma distribuição de probabilidades que considera graus de liberdade (\(\nu\), letra grega ni ).
Sob \(H_0\) é semelhante à distribuição normal padronizada, centrada em \(t=0\), mas com suas caudas mais “pesadas” (desvio-padrão > 1). Não é uma única curva, mas uma família delas, variando os graus de liberdade: podemos pensar na distribuição \(t\) como um avanço histórico em relação à distribuição normal padronizada - em um teste \(z\) o desvio-padrão populacional é conhecido; no teste \(t\) usa-se o desvio-padrão amostral como estimador. A incerteza adicional pela falta de conhecimento do desvio-padrão populacional é considerada através dos graus de liberdade, alterando a distribuição sobre a qual a estatística do teste funciona.
Os graus de liberdade dependem do tamanho da amostra; quanto menor, mais pesadas são as caudas e, portanto, diferenças numéricas precisam que ser maiores para que consigamos rejeitar a hipótese nula.
|
Experimente Animacao_t_central.R para ver o aspecto da distribuição \(t\) e observe:
|
R dispõe de uma pequena família de funções básicas para cada tipo de distribuição. Estes pequenos conjuntos são análogos uns aos outros conjuntos, facilitando o aprendizado.
A família de funções para a distribuição binomial é:
que dependem do número de eventos (size) e sucesso de cada evento (prob). Esta distribuição serve, por exemplo, para modelar uma moeda. Por exemplo, uma moeda bem balanceada tem \(50\%\) de probabilidade de sair coroa. Caso fosse testada com 6 jogadas, esperamos que 3 coroas sejam o mais provável, e o gráfico com a distribuição de probabilidades pode usar o seguinte código:
jogadas <- seq(0:6)
probabilidades <- dbinom(x=jogadas, size=6, prob=0.5)
plot(jogadas, probabilidades,
xlab="coroas", ylab="probabilidade",
xlim=c(0,6), type="h")Para uma moeda desbalanceada, viciada para dar \(2 \over 3\) de coroas, esperamos 4 coroas em 6 jogadas, como vemos alterando o parâmetro prob:
jogadas <- seq(0:6)
probabilidades <- dbinom(x=jogadas, size=6, prob=2/3)
plot(jogadas, probabilidades,
xlab="coroas", ylab="probabilidade",
xlim=c(0,6), type="h", lty=2)A curva produzida para a moeda desbalanceada é a mesma daquela obtida pela moeda balanceada, apenas transladada para a direita. A dificuldade para distinguir uma moeda da outra está em observar que a moeda de \(50\%\) tem certa probabilidade de sair com 4 coroas, bem como a desbalanceada pode obter 3.
A distribuição normal tem conjunto de funções similar às da binomial:
na qual necessitamos da média (mean) e desvio-padrão (sd) para sua caracterização. É usada com variáveis quantitativas contínuas, portanto os gráficos devem usar linhas contínuas.
Por comparação exploraremos distribuições normais com médias de 3 e 4 e desvio-padrão de 1.22.
|
Este desvio-padrão não foi escolhido ao acaso. O desvio-padrão de uma binomial é dado por: \(sd_{binomial} = \sqrt{ n \cdot p \cdot (1-p) }\) onde \(n\) é o número de jogadas, \(p\) é a probabilidade do evento. Então, para o exemplo acima, \(sd_{binomial} = \sqrt{ 6 \cdot 0.5 \cdot (1-0.5) } \approx 1.22\). Para este exemplo procuramos mostrar uma distribuição normal com formato similar à binomial do exemplo anterior. |
O código para ilustrar distribuições normais é muito similar aos mostrados para as distribuições binomais, mas usando a função dnorm() no lugar de dbinom().
Para mostrar uma distribuição normal com média de 3 podemos usar:
valores <- seq(from=0, to=6, by=0.01)
densidades <- dnorm(x=valores, mean=3, sd=1.22)
plot(valores, densidades,
xlab="valor", ylab="densidade",
type="l")e para média de 4:
valores <- seq(from=0, to=6, by=0.01)
densidades <- dnorm(x=valores, mean=4, sd=1.22)
plot(valores, densidades,
xlab="valor", ylab="densidade",
type="l", lty=2)Observamos, novamente, a curva transladada para a direita quando a média aumenta de 3 para 4. O problema estatístico é análogo a distinguir uma moeda balanceada de uma viciada em \(2 \over 3\): saber se duas distribuições com médias numericamente diferentes podem ser tratadas como diversas.
A decisão estatística, considerando \(H_0: \mu_A = \mu_B\), é tomada com base nas distribuições normais padronizadas: é o caso de um teste \(z\), e as duas distribuições exemplificadas aqui corresponderiam aproximadamente a:
z <- seq(from=-4*1.22, to=4*1.22, by=0.01)
H0_z <- dnorm(x=z, mean=0, sd=1)
plot(z, H0_z,
xlab="z", ylab="densidade",
type="l")
delta <- (4-3)/1.22
H1_z <- dnorm(x=z, mean=delta, sd=1)
lines(z, H1_z, lty=2)Esta representação ilustra normais padronizadas no intervalo de \(\pm 4~\text{desvios-padrão}\). A distribuição de referência, que tinha média de 4 (correspondendo a \(H_0: \mu_A = \mu_B\), linha sólida) foi centrada em zero. A distribuição correspondente à média de 4 (\(H_1: \mu_A \ne \mu_B\), linha pontilhada) está a aproximadamente a 0.82 unidades de desvio-padrão acima da média de referência (\((4-3)/1.22 \approx 0.8196\)).
Para a distribuição \(t\) as funções são:
Como as distribuições normais, as funções t são para variáveis quantitativas contínuas. Já são padronizadas e, portanto, não aparece média e desvio-padrão entre seus parâmetros. Em vez disto, aparecem duas novidades:
No caso, df é relacionado com o tamanho da amostra (\(df = n-1\)), e ncp define a translação da distribuição \(t\), representando \(H_1\). A curva correspondente a \(H_0\) tem ncp = 0.
O código similar ao das normais padronizadas é:
t <- seq(from=-6, to=6, by=0.01)
H0_t <- dt(x=t, df=10-1, ncp=0)
plot(t, H0_t,
xlab="t", ylab="densidade",
type="l")
H1_t <- dt(x=t, df=10-1, ncp=0.82)
lines(t, H1_t, lty=2)Algo além da translação ocorreu quando ncp não é zero: observe que a curva pontilhada é ligeiramente mais baixa que a de linha sólida. Mais difícil de perceber é que a curva pontilhada é assimétrica. É mais visível com valor maior de ncp:
t <- seq(from=-6, to=6, by=0.01)
H0_t <- dt(x=t, df=10-1, ncp=0)
plot(t, H0_t,
xlab="t", ylab="densidade",
type="l")
H1_t <- dt(x=t, df=10-1, ncp=2)
lines(t, H1_t, lty=2)Na estatística \(t\), então, as curvas que representam a hipótese alternativa são transladas mas, também, assimétricas.
| O código disponível em Animacao_t_nao_central.R compara a curva observada sob \(H_0\) na animação anterior com diversas curvas representando \(H_1\)s. Observe, sob \(H_1\), assimetria das distribuições, e as áreas correspondentes a \(\beta\) e ao poder do teste (\(1 - \beta\)). |
Análise estatística inferencial é o processo de estimar características de uma população a partir de uma amostra, através de teste da hipótese nula, usando seus estimadores.
Suspeita-se de que um medicamento vasodilatador (Nifedipina) para Hipertensão Arterial, amplamente receitado, esteja aumentando a frequência cardíaca dos pacientes.
É sabido que a frequência cardíaca na população normal tem Distribuição Normal com média 70 bpm.
Para verificar essa suspeita, planejou-se obter uma amostra aleatória de 50 pacientes que recebem Nifedipina para se medir a frequência cardíaca.
\(H_0: \mu_{nifedipina} = \mu_0\)
\(H_1: \mu_{nifedipina} > \mu_0\)
Adota-se \(\mu_0 = 70~\text{bpm}\)
|
O teste é unicaudal (só investigamos se há aumento da frequência cardíaca) e a direção é explícita em \(H_1\), com o símbolo \(>\). É comum encontrar a anotação da hipótese nula com o símbolo complementar (\(\le\) neste exemplo): \(H_0: \mu_{nifedipina} \le \mu_0\) \(H_1: \mu_{nifedipina} > \mu_0\) Mas optamos por usar o símbolo de igualdade (\(=\)) porque mais adequadamente espelha o que se espera de \(H_0\), a ausência de efeito. Matematicamente, também, é equivalente (GATÁS RR (1978, p. 220-223) Elementos de Probabilidade e Inferência. SP: Atlas.) |
A amostra de 50 pacientes forneceu:
72, 74, 70, 70, 69, 71, 72, 71, 69, 74, 71, 71, 70, 73, 69, 68, 68, 71, 71, 72, 70, 69, 73, 69, 71, 70, 72, 73, 70, 72, 67, 72, 67, 68, 69, 72, 70, 70, 70, 71, 74, 67, 69, 71, 71, 73, 71, 71, 70, 71
bpm <- c(72, 74, 70, 70, 69, 71, 72, 71, 69, 74, 71, 71, 70, 73, 69, 68, 68,
71, 71, 72, 70, 69, 73, 69, 71, 70, 72, 73, 70, 72, 67, 72, 67, 68,
69, 72, 70, 70, 70, 71, 74, 67, 69, 71, 71, 73, 71, 71, 70, 71)
sumario <- summary(bpm)
print(sumario) Min. 1st Qu. Median Mean 3rd Qu. Max.
67.00 69.25 71.00 70.58 72.00 74.00
par(mfrow=c(2,2))
hist(bpm, freq=TRUE,
main="Dados amostrais",
xlab="Batimentos cardíacos", ylab="Contagem")
divisoes <- seq(from=min(bpm),to=max(bpm),by=(max(bpm)-min(bpm))/5)
hist(bpm, freq=TRUE, breaks=divisoes,
main="Dados amostrais",
xlab="Batimentos cardíacos", ylab="Contagem")
divisoes <- seq(from=min(bpm),to=max(bpm),by=(max(bpm)-min(bpm))/8)
hist(bpm, freq=TRUE, breaks=divisoes,
main="Dados amostrais",
xlab="Batimentos cardíacos", ylab="Contagem")
divisoes <- seq(from=min(bpm),to=max(bpm),by=(max(bpm)-min(bpm))/9)
hist(bpm, freq=TRUE, breaks=divisoes,
main="Dados amostrais",
xlab="Batimentos cardíacos", ylab="Contagem")par(mfrow=c(1,1))par(mfrow=c(1,2))
boxplot (bpm,
main="Dados amostrais",
ylab="Batimentos cardíacos", xlab="")
boxplot (bpm, horizontal = TRUE,
main="Dados amostrais",
xlab="Batimentos cardíacos", ylab="")par(mfrow=c(1,1))densprob <- density(bpm)
plot (densprob,
main="Distribuição dos dados amostrais",
xlab="Batimentos cardíacos", ylab="Densidade")no “olhômetro”, parece aproximar-se da distribuição normal?
# dados
bpm <- c(72, 74, 70, 70, 69, 71, 72, 71, 69, 74, 71, 71, 70, 73, 69, 68, 68,
71, 71, 72, 70, 69, 73, 69, 71, 70, 72, 73, 70, 72, 67, 72, 67, 68,
69, 72, 70, 70, 70, 71, 74, 67, 69, 71, 71, 73, 71, 71, 70, 71)
# density plot
densprob <- density(bpm)
plot (densprob,
main="Distribuição dos dados amostrais",
xlab="Batimentos cardíacos", ylab="Densidade")
# distribuicao com media +- 4 desvios-padrao
media_bpm <- mean(bpm, na.rm = TRUE)
dp_bpm <- sd(bpm, na.rm = TRUE)
x <- seq(media_bpm-4*dp_bpm,media_bpm+4*dp_bpm,by=0.1)
distnormal <- dnorm(x, mean=media_bpm, sd=dp_bpm)
lines(x, distnormal, lty=2)Para um teste \(t\) para uma condição, unilateral à direita (note o uso de “greater”), é necessário fornecer o valor de referência (mu_pop <- 70) executado com:
bpm <- c(72, 74, 70, 70, 69, 71, 72, 71, 69, 74, 71, 71, 70, 73, 69, 68, 68,
71, 71, 72, 70, 69, 73, 69, 71, 70, 72, 73, 70, 72, 67, 72, 67, 68,
69, 72, 70, 70, 70, 71, 74, 67, 69, 71, 71, 73, 71, 71, 70, 71)
mu_pop <- 70
alfa <- 0.05
t_out <- t.test(bpm, mu=mu_pop,
conf.level = 1-alfa, alternative = "greater")
print (t_out)
One Sample t-test
data: bpm
t = 2.312, df = 49, p-value = 0.01251
alternative hypothesis: true mean is greater than 70
95 percent confidence interval:
70.15942 Inf
sample estimates:
mean of x
70.58
Para a decisão estatística, observe o valor da estatística do teste, guardada em t_out$statistic=2.3120489 e o valor-\(p\) associado em t_out$p.value=0.0125097.
Para \(\alpha=0.05\), rejeita-se \(H_0\) e, portanto, o uso de nifedipina está associada ao aumento da frequência cardíaca neste estudo.
Para \(\alpha=0.01\), não se rejeita \(H_0\) e, portanto, não há elementos neste estudo para afirmar-se que o uso de nifedipina está associada ao aumento da frequência cardíaca.
FALTOU escolher \(\alpha\) no planejamento do estudo.
Disponibilizamos o RScript Nifedipina.R que reúne os procedimentos acima, gerando alguns resultados e gráficos adicionais.
BPM :
Min. 1st Qu. Median Mean 3rd Qu. Max.
67.00 69.25 71.00 70.58 72.00 74.00
Desvio-padrao = 1.77
Tamanho da amostra = 50
One Sample t-test
data: Nifedipina$BPM
t = 2.312, df = 49, p-value = 0.01251
alternative hypothesis: true mean is greater than 70
95 percent confidence interval:
70.15942 Inf
sample estimates:
mean of x
70.58
Referencial teorico:
valor-p = 0.01250973
d de Cohen = 0.3269731 (Pequeno)
sob H0: distribuicao t com media = 0 e d.p. = 1.021055
sob H1: distribuicao t com media = 2.348206 e d.p. = 1.049534
Note que:
|
A outra maneira de decidir é utilizar o intervalo de confiança (IC). Note que a função t.test() exigiu o parâmetro \(\alpha\); necessário, justamente, para o cálculo do IC. No exemplo, fornecemos \(\alpha=0.05\) e foi computado IC95 = [70.1594209, Inf] (guardado na variável t_out$conf.int). O limite superior é infinito (Inf significa \(\infty\)) porque o teste é unilateral. O valor populacional (\(\mu=70~bpm\)) está fora do intervalo, levando à rejeição de \(H_0\) para este alfa. Caso o teste fosse executado com \(\alpha=0.01\) teríamos:
O valor populacional (\(\mu=70~bpm\)) agora está dentro do intervalo, levando à não rejeição de \(H_0\).
|
Alterando-se, no início de Nifedipina.R o valor de alfa:
alfa <- 0.01 # nivel de significancia adotadoa saída altera-se de acordo, mostrando \(p>\alpha\) e a não-rejeição de \(H_0\):
BPM :
Min. 1st Qu. Median Mean 3rd Qu. Max.
67.00 69.25 71.00 70.58 72.00 74.00
Desvio-padrao = 1.77
Tamanho da amostra = 50
One Sample t-test
data: Nifedipina$BPM
t = 2.312, df = 49, p-value = 0.01251
alternative hypothesis: true mean is greater than 70
99 percent confidence interval:
69.97671 Inf
sample estimates:
mean of x
70.58
Referencial teorico:
valor-p = 0.01250973
d de Cohen = 0.3269731 (Pequeno)
sob H0: distribuicao t com media = 0 e d.p. = 1.021055
sob H1: distribuicao t com media = 2.348206 e d.p. = 1.049534
\(H_0\) é a hipótese nula (sempre abrange a igualdade). \(H_1\) é a hipótese alternativa (oposta da \(H_0\) ). \(H_0\) pode ser não rejeitada ou rejeitada (a rejeição deve ser baseada em evidências obtidas a partir da amostra).
A decisão estatística sempre envolve possíveis erros: que são:
A decisão sobre o teste depende da comparação entre a probabilidade de se observar uma diferença sob a hipótese nula, dada uma amostra de tamanho \(n\) e a probabilidade do erro do tipo I (\(\alpha\)) escolhida previamente:
Também chamado de teste \(t\) relacionado, aplica-se tipicamente às situações em que o mesmo indivíduo (ou a mesma unidade experimental) tem uma variável quantitativa medida em dois momentos ou duas condições experimentais.
A pesquisadora Yob está interessada na violência de massa durante as partidas de futebol. Ela pensa que a violência do grupo é resultado dos assentos desconfortáveis do estádio.
Por isso, Yob modifica dois estádios diferentes na Inglaterra. Em um estádio coloca assentos bem apertados e desconfortáveis. No outro, instala assentos confortáveis, com muito espaço para as pernas e entre os assentos adjacentes.
A professora organiza uma competição, de modo que um clube jogue metade das partidas em um estádio e a outra metade no outro estádio. Ela prevê que o número de prisões e expulsões será maior no estádio que apresenta os assentos mais desconfortáveis.
Ela acompanha um grupo de 12 fãs adolescentes agressivos e grosseiros do clube e registra o número de vezes que cada um é preso ou expulso do estádio.
Aqui disponibilizamos o RScript Violencia_estadios.R e destacamos seus principais trechos.
Os dados estão disponíveis na planilha Excel Violencia_estadios.xlsx:
library(readxl)
Dtfrm <- read_excel("Violencia_estadios.xlsx", sheet = "dependente")
# diferenca entre as condições experimentais
Dtfrm$dif <- Dtfrm$Conforto - Dtfrm$Desconforto
print(Dtfrm)# A tibble: 12 x 4
Adolescente Desconforto Conforto dif
<chr> <dbl> <dbl> <dbl>
1 a 8 3 -5
2 b 5 2 -3
3 c 4 4 0
4 d 6 6 0
5 e 4 2 -2
6 f 8 1 -7
7 g 9 6 -3
8 h 10 3 -7
9 i 7 4 -3
10 j 8 1 -7
11 k 6 4 -2
12 l 7 3 -4
Abra a planilha para verificar como os dados foram guardados e como aparecem ao serem lidos. Note que read_excel() lê a aba “dependente” da planilha, na qual cada linha tem as observações de um indivíduo, feitas nas duas condições experimentais.
Exibe a estatística descritiva:
cat("\nEsquema de 5 estatísticas de Tukey & média\n")
cat("\nTamanho da amostra: ",length(Dtfrm$dif),"\n", sep="")
cat("\nNúmero de ocorrências em ",names(Dtfrm)[2],":\n", sep="")
sumario <- summary(Dtfrm[[2]], digits = 3)
print (sumario)
cat("\nNúmero de ocorrências ",names(Dtfrm)[3],":\n", sep="")
sumario <- summary(Dtfrm[[3]], digits = 3)
print (sumario)
cat("Diferença do número de ocorrências (",names(Dtfrm)[3]," - ",names(Dtfrm)[2],"):\n", sep="")
sumario <- summary(Dtfrm$dif, digits = 3)
print (sumario)
Esquema de 5 estatísticas de Tukey & média
Tamanho da amostra: 12
Número de ocorrências em Desconforto:
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.00 5.75 7.00 6.83 8.00 10.00
Número de ocorrências Conforto:
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 2.00 3.00 3.25 4.00 6.00
Diferença do número de ocorrências (Conforto - Desconforto):
Min. 1st Qu. Median Mean 3rd Qu. Max.
-7.00 -5.50 -3.00 -3.58 -2.00 0.00
Gera alguns gráficos para estatística descritiva (abra o código R para ver como os gráficos foram gerados):
Faz testes de estatística inferencial:
corr_test <- cor.test(Dtfrm$Desconforto,Dtfrm$Conforto)
cat("\nCoeficiente de correlacao de Pearson\n")
print (corr_test)
Coeficiente de correlacao de Pearson
Pearson's product-moment correlation
data: Dtfrm$Desconforto and Dtfrm$Conforto
t = 0.04565, df = 10, p-value = 0.9645
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5641406 0.5835022
sample estimates:
cor
0.01443426
mostrando que os comportamentos do mesmo indivíduo em cada uma das condições experimentais não estão associados.
cat("Analise de significancia estatistica: valor-p\n")
t_out <- t.test(Dtfrm$dif,mu=0,alternative="less")
print(t_out)Analise de significancia estatistica: valor-p
One Sample t-test
data: Dtfrm$dif
t = -4.9592, df = 11, p-value = 0.0002147
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
-Inf -2.285695
sample estimates:
mean of x
-3.583333
mostrando que há diferença entre as condições experimentais (rejeição de \(H_0\), valor-\(p < \alpha\), intervalo de confiança 95% não inclui o valor zero). O teste foi feito com \(\text{Conforto} -\text{Desconforto}\), encontrando-se números negativos: então NPE é maior em \(\text{Desconforto}\); o número de NPE é significantemente maior nos estádios desconfortáveis.
Computa, também, valores para a significância prática:
F <- t^2
df <- t_out$parameter
eta2 <- F/(F+df)
# Elis P (2010) The essential guide to effect sizes. Cambrige
if (eta2 <0.01) {mag_eta2<-c("Desprezivel")}
if (eta2>=0.01 && eta2<0.06) {mag_eta2<-c("Pequeno")}
if (eta2>=0.06 && eta2<0.14) {mag_eta2<-c("Intermediario")}
if (eta2>=0.14) {mag_eta2<-c("Grande")}
R2aj <- (F-1)/(F+df)
# tamanho de efeito d de Cohen
dp <- sd(Dtfrm$dif)
m <- t_out$estimate
d <- abs(t_out$statistic)/sqrt(t_out$parameter+1)
# Sawilowsky, S (2009) New effect size rules of thumb. Journal of Modern Applied Statistical Methods 8(2): 467-74.
if (d<0.01) {mag_Cohen<-c("Desprezivel")}
if (d>=0.01 && d<0.2) {mag_Cohen<-c("Muito pequeno")}
if (d>=0.2 && d<0.5) {mag_Cohen<-c("Pequeno")}
if (d>=0.5 && d<0.8) {mag_Cohen<-c("Intermediario")}
if (d>=0.8 && d<1.2) {mag_Cohen<-c("Grande")}
if (d>=1.2 && d<2) {mag_Cohen<-c("Muito grande")}
if (d>=2) {mag_Cohen<-c("Enorme")}
cat("Analise de significancia pratica: tamanho de efeito\n")
cat("\td de Cohen = ",d," (",mag_Cohen,")","\n",sep="")
cat("\teta^2 = R^2 = ",eta2," (",mag_eta2,")\n",sep="")Analise de significancia pratica: tamanho de efeito
d de Cohen = 1.431599 (Muito grande)
eta^2 = R^2 = 0.3270348 (Grande)
A tabela implementada no RScript para a intensidade do efeito calculado por eta^2 (\(\eta^2\)) é: O SNAP-Ed (Supplemental Nutrition Assistance Program Education) é um programa baseado em evidências que ajuda as pessoas a terem uma vida mais saudável.
O SNAP-Ed ensina às pessoas que usam ou qualificam para o SNAP uma boa nutrição e como fazer com que o seu dinheiro de alimentação se estenda ainda mais.
Os participantes do SNAP-Ed também aprendem a ser fisicamente ativos.
Brendon Small e Coach McGuirk fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as duas turmas.
Desenvolvemos Nutricao.R para as análises descritiva e inferencial. Abaixo destacamos seus techos principais. Neste RScript aproveitamos funções de alguns pacotes (verifique se estão instalados em seu computador)
library(readxl)
library(car)Loading required package: carData
library(lattice)
library(ggplot2)Registered S3 methods overwritten by 'ggplot2':
method from
[.quosures rlang
c.quosures rlang
print.quosures rlang
library(rcompanion)Os dados estão na planilha Nutricao.xlsx:
Dtfrm <- read_excel("Nutricao.xlsx")
# os instrutores devem ser tratados como fator
Dtfrm$Instructor <- as.factor(Dtfrm$Instructor)
Dtfrm$Instructor <- factor(Dtfrm$Instructor, levels=unique(Dtfrm$Instructor))
print(Dtfrm)# A tibble: 40 x 3
Instructor Student Sodium
<fct> <chr> <dbl>
1 Brendon Small a 1200
2 Brendon Small b 1400
3 Brendon Small c 1350
4 Brendon Small d 950
5 Brendon Small e 1400
6 Brendon Small f 1150
7 Brendon Small g 1300
8 Brendon Small h 1325
9 Brendon Small i 1425
10 Brendon Small j 1500
# … with 30 more rows
Verificamos se os dados estão coerentes:
res_sodiumbyinstructor <- summary.data.frame(Dtfrm, digits=2)
print (res_sodiumbyinstructor) Instructor Student Sodium
Brendon Small:20 Length:40 Min. : 950
Coach McGuirk:20 Class :character 1st Qu.:1150
Mode :character Median :1250
Mean :1267
3rd Qu.:1362
Max. :1700
Gráficos sugeridos:
grf <- boxplot(Sodium ~ Instructor, data=Dtfrm, ylab=names(Dtfrm)[3],
xlab="Instructor")print(grf)$stats
[,1] [,2]
[1,] 950 1000.0
[2,] 1150 1137.5
[3,] 1300 1212.5
[4,] 1400 1350.0
[5,] 1700 1525.0
$n
[1] 20 20
$conf
[,1] [,2]
[1,] 1211.675 1137.424
[2,] 1388.325 1287.576
$out
numeric(0)
$group
numeric(0)
$names
[1] "Brendon Small" "Coach McGuirk"
grf <- lattice::xyplot(Sodium ~ Instructor, data=Dtfrm, type=c("p","a"), jitter.x=TRUE)
print(grf)grf <- car::densityPlot(Sodium~Instructor, data=Dtfrm, rug=TRUE)print(grf)$`Brendon Small`
Call:
adaptiveKernel(x = x[g == group], bw = if (is.numeric(bw)) bw[group] else bw, adjust = adjust[group])
Data: x[g == group] (500 obs.); Bandwidth 'bw' = 92.23
x y
Min. : 673.3 Min. :3.574e-05
1st Qu.: 999.2 1st Qu.:1.735e-04
Median :1325.0 Median :4.381e-04
Mean :1325.0 Mean :7.619e-04
3rd Qu.:1650.8 3rd Qu.:1.218e-03
Max. :1976.7 Max. :2.451e-03
$`Coach McGuirk`
Call:
adaptiveKernel(x = x[g == group], bw = if (is.numeric(bw)) bw[group] else bw, adjust = adjust[group])
Data: x[g == group] (500 obs.); Bandwidth 'bw' = 70.4
x y
Min. : 788.8 Min. :2.804e-05
1st Qu.:1025.6 1st Qu.:2.137e-04
Median :1262.5 Median :7.231e-04
Mean :1262.5 Mean :1.050e-03
3rd Qu.:1499.4 3rd Qu.:1.708e-03
Max. :1736.2 Max. :3.032e-03
SumMM <- rcompanion::groupwiseMean(Sodium ~ Instructor,
data = Dtfrm,
conf = 0.95,
digits = 3,
traditional = FALSE,
percentile = TRUE)
grf <- ggplot2::ggplot(SumMM, ggplot2::aes(x = Instructor, y = Mean)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin = Percentile.lower,
ymax = Percentile.upper),
width = 0.05, size = 0.5) +
ggplot2::geom_point(shape = 15,
size = 4) +
ggplot2::theme_bw() +
ggplot2::theme(axis.title = ggplot2::element_text(face = "bold")) +
ylab(names(Dtfrm)[3])
print(grf)|
Repare o que acontece com a linha
|
Definimos alfa:
alfa <- 0.05 # nivel de significancia adotadoO teste \(t\) a ser aplicado é de Satterthwaite (apesar do R exibir como teste de Welch), conforme as seguintes referências:
Para operacionalizar o teste, calculamos:
# separa os dois instrutores
SodiumBS <- subset(Dtfrm,select=Sodium,subset=Instructor=="Brendon Small",drop=TRUE)
SodiumCM <- subset(Dtfrm,select=Sodium,subset=Instructor=="Coach McGuirk",drop=TRUE)
# dados da amostra
nA <- sum(!is.na(SodiumBS))
nB <- sum(!is.na(SodiumCM))
# significancia estatistica
t_out <- t.test(Sodium ~ Instructor, data = Dtfrm)
# significancia pratica
t <- t_out$statistic # estatistica de teste t
df <- t_out$parameter # graus de liberdade
dfefic <- (df-min(nA,nB))/(nA+nB-2-min(nA,nB))exibimos o resultado:
cat("\nTamanho das amostras: \n", sep="")
cat ("\tBrendon Small: n = ", nA, "\n", sep="")
cat ("\tCoach McGuirk: n = ", nB, "\n", sep="")
cat ("\n")
cat("Analise de significancia estatistica: valor-p\n")
cat("Teste t de Satterthwaite\n")
print(t_out)
cat("Eficiencia do numero de graus de liberdade =",dfefic,"\n\n")
Tamanho das amostras:
Brendon Small: n = 20
Coach McGuirk: n = 20
Analise de significancia estatistica: valor-p
Teste t de Satterthwaite
Welch Two Sample t-test
data: Sodium by Instructor
t = 0.76722, df = 34.893, p-value = 0.4481
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-67.91132 150.41132
sample estimates:
mean in group Brendon Small mean in group Coach McGuirk
1287.50 1246.25
Eficiencia do numero de graus de liberdade = 0.8273926
e construímos os gráficos:
# distribuicao t sob H0 (central: ncp = 0)
tH0 <- rt(1e6, df)
dtH0 <- density(tH0)
# distribuicao t sob H1, ncp = t
# ncp
F <- t^2 # estatistica de teste F de Fisher
eta2 <- F/(F+df)
f2 <- eta2/(1-eta2) # f de Cohen
ncp <- df*f2 # parametro de nao-centralidade: ncp = F
tH1 <- rt(1e6, df, ncp)
dtH1 <- density(tH1)
# media e dp das dist. t central e nao-central
mediaH0 <- 0
dpH0 <- sqrt(df/(df-2))
beta <- sqrt(df/2)*gamma((df-1)/2)/gamma(df/2)
mediaH1 <- ncp*beta
dpH1 <- sqrt((df*(1+ncp^2)/(df-2))-mediaH1^2)cat("\nReferencial teorico:\n")
cat("\tvalor-p =",t_out$p.value,"\n")
cat("sob H0: distribuicao t com media = ",mediaH0," e d.p. = ",dpH0,"\n",sep="")
cat("sob H1: distribuicao t com media = ",mediaH1," e d.p. = ",dpH1,"\n",sep="")
cat("\n")
Referencial teorico:
valor-p = 0.4481092
sob H0: distribuicao t com media = 0 e d.p. = 1.029953
sob H1: distribuicao t com media = 0.6016735 e d.p. = 1.032641
e construímos os gráficos:
# graficos
for (g in 1:2)
{
# limites de x
min_x <- min(mediaH0-3*dpH0, mediaH1-3*dpH1)
max_x <- max(mediaH0+3*dpH0, mediaH1+3*dpH1)
if (g == 1)
{
plot(dtH0,
main=paste("Teste t independente\ndf =",round(df,3),", t =",round(t,5),", alfa =",alfa),
xlab="t",
xlim=c(min_x,max_x),
lwd=1, lty=1
)
}
if (g == 2)
{
plot(dtH0,
main=NA,
xlab="t",
xlim=c(min_x,max_x),
lwd=1, lty=1
)
}
qalfa <- qt(c(alfa/2,1-alfa/2),df,0)
abline(v=qalfa[1], lty = 3)
abline(v=qalfa[2], lty = 3)
if (g==1)
{
abline(v=abs(t),lwd=1,lty=2)
abline(v=-abs(t),lwd=1,lty=2)
# area do valor p
polx <- dtH0$x[dtH0$x>=abs(t)]; polx <- c(min(polx),polx,max(polx))
poly <- dtH0$y[dtH0$x>=abs(t)]; poly <- c(0,poly,0)
polygon(polx,poly,border="#EE802622",col="#EE802688",lwd=5)
polx <- dtH0$x[dtH0$x<=-abs(t)]; polx <- c(min(polx),polx,max(polx))
poly <- dtH0$y[dtH0$x<=-abs(t)]; poly <- c(0,poly,0)
polygon(polx,poly,border="#EE802622",col="#EE802688",lwd=5)
}
if (g==2)
{
# H1
lines(dtH1,lwd=3,lty=1)
# area alfa
polx <- dtH0$x[dtH0$x<=qalfa[1]]; polx <- c(min(polx),polx,max(polx))
poly <- dtH0$y[dtH0$x<=qalfa[1]]; poly <- c(0,poly,0)
polygon(polx,poly,border="#a3261b22",col="#a3261b88",lwd=5)
polx <- dtH0$x[dtH0$x>=qalfa[2]]; polx <- c(min(polx),polx,max(polx))
poly <- dtH0$y[dtH0$x>=qalfa[2]]; poly <- c(0,poly,0)
polygon(polx,poly,border="#a3261b22",col="#a3261b88",lwd=5)
# area beta
polx <- dtH1$x[dtH1$x>=qalfa[1] & dtH1$x<=qalfa[2]]; polx <- c(min(polx),polx,max(polx))
poly <- dtH1$y[dtH1$x>=qalfa[1] & dtH1$x<=qalfa[2]]; poly <- c(0,poly,0)
polygon(polx,poly,border="#4EB26522",col="#4EB26588",lwd=8)
}
# legenda
if (g==1)
{
p_txt <- t_out$p.value;
if (p_txt > 0.001)
{
p_txt <- round(p_txt,4)
} else
{
p_txt <- format(format(p_txt, scientific = TRUE, digits = 4))
}
legend ("topright",
c("H0","t obs.","t crit.",paste("p=",p_txt,sep="")),
lwd=c(1,1,1,10),
lty=c(1,2,3,1),
pch=NA,
col=c("black","black","black","#EE802688"),
box.lwd=0, bg="transparent")
}
if (g==2)
{
legend ("topright",
c("H0","H1","t crit.","alfa","beta"),
lwd=c(1,3,1,10,10),
lty=c(1,1,3,1,1),
pch=NA,
col=c("black","black","black","#a3261b88","#4EB26588"),
box.lwd=0, bg="transparent")
}
}Verifique e interprete a saída. Não rejeitamos \(H_0\): não há elementos para afirmar que há diferença de resultado, quanto à ingestão de sódio, quando comparamos os dois grupos submetidos a diferentes programas educacionais.
|
Existe um outro teste \(t\) feito por bootstrapping (do package MKinfer):
A novidade são as primeiras linhas. O que aparece após “Results without bootstrap” é o teste de Welch visto acima.
|
Calculamos os tamanhos de efeito:
# Elis P (2010) The essential guide to effect sizes. Cambrige
if (eta2 <0.01) {mag_eta2<-c("Desprezivel")}
if (eta2>=0.01 && eta2<0.06) {mag_eta2<-c("Pequeno")}
if (eta2>=0.06 && eta2<0.14) {mag_eta2<-c("Intermediario")}
if (eta2>=0.14) {mag_eta2<-c("Grande")}
d <- abs(t)/sqrt(1/((1/nA)+(1/nB)))
# Sawilowsky, S (2009) New effect size rules of thumb. Journal of Modern Applied Statistical Methods 8(2): 467-74.
if (d<0.01) {mag_Cohen<-c("Desprezivel")}
if (d>=0.01 && d<0.2) {mag_Cohen<-c("Muito pequeno")}
if (d>=0.2 && d<0.5) {mag_Cohen<-c("Pequeno")}
if (d>=0.5 && d<0.8) {mag_Cohen<-c("Intermediario")}
if (d>=0.8 && d<1.2) {mag_Cohen<-c("Grande")}
if (d>=1.2 && d<2) {mag_Cohen<-c("Muito grande")}
if (d>=2) {mag_Cohen<-c("Enorme")}
g <- d*(1-3/(4*df-1))
if (g<0.01) {mag_Hedges<-c("Desprezivel")}
if (g>=0.01 && g<0.2) {mag_Hedges<-c("Muito pequeno")}
if (g>=0.2 && g<0.5) {mag_Hedges<-c("Pequeno")}
if (g>=0.5 && g<0.8) {mag_Hedges<-c("Intermediario")}
if (g>=0.8 && g<1.2) {mag_Hedges<-c("Grande")}
if (g>=1.2 && g<2) {mag_Hedges<-c("Muito grande")}
if (g>=2) {mag_Hedges<-c("Enorme")}
# selecao de modelo
R2aj <- (F-1)/((F-1)+df+1)
omega2 <- (F-1)/((F-1)+df+2)e exibimos o resultado:
cat("Analise de significancia pratica: tamanho de efeito\n")
cat("\td de Cohen = ",d," (",mag_Cohen,")","\n",sep="")
cat("\tg de Hedges = ",g," (",mag_Hedges,")","\n",sep="")
cat("\teta^2 = R^2 = ",eta2," (",mag_eta2,")\n",sep="")
cat("\nSelecao de modelo:\n")
cat("\tR^2 ajustado =",R2aj,"\n")
cat("\tomega^2 = ", omega2,"\n")Analise de significancia pratica: tamanho de efeito
d de Cohen = 0.2426174 (Pequeno)
g de Hedges = 0.2373649 (Pequeno)
eta^2 = R^2 = 0.01658973 (Pequeno)
Selecao de modelo:
R^2 ajustado = -0.01159381
omega^2 = -0.01127601
É muito comum, em publicações, que somente tenhamos acesso às medidas-resumo (número de participantes, média, desvio-padrão e correlação). Nestes casos, os RScripts acima não são utilizáveis.
Para fazer os testes \(t\) de Welch ( independentes) ou testes \(t\) relacionados (pareados), quando os dados brutos não estão disponíveis, criamos os seguintes scripts:
Estude e aprenda a modificar estes RScripts para seu uso.
Um código que incorpora o teste \(t\) e um cálculo de d de Cohen (requer o package lsr) mais adequado à versão robusta do teste \(t\) de Welch é:
library(readxl)
library(MKinfer)
library(lsr)
Dados <- readxl::read_excel("Nutricao.xlsx")
MKinfer::boot.t.test(Sodium ~ Instructor, data=Dados, R=10^6)
Bootstrapped Welch Two Sample t-test
data: Sodium by Instructor
bootstrapped p-value = 0.447
95 percent bootstrap percentile confidence interval:
-61.25 143.75
Results without bootstrap:
t = 0.76722, df = 34.893, p-value = 0.4481
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-67.91132 150.41132
sample estimates:
mean in group Brendon Small mean in group Coach McGuirk
1287.50 1246.25
# d de Cohen para o teste t de Welch
d <- lsr::cohensD(Sodium ~ Instructor, data=Dados, method="unequal")
# Sawilowsky, S (2009) New effect size rules of thumb.
# Journal of Modern Applied Statistical Methods, 8(2): 467-74.
if (0 <= d && d < 0.01) {dc <- "negligible"}
if (0.01 <= d && d < 0.2) {dc <- "very small"}
if (0.2 <= d && d < 0.5) {dc <- "small"}
if (0.5 <= d && d < 0.8) {dc <- "medium"}
if (0.8 <= d && d < 1.2) {dc <- "large"}
if (1.2 <= d && d < 2) {dc <- "very large"}
if (2 <= d && d < Inf) {dc <- "huge"}
cat("Cohen's d = ", d, " (",dc, " effect size)", sep="")Cohen's d = 0.2426174 (small effect size)
Um exemplo caricato (apenas 4 medidas em cada grupo) é dado por:
estatm <- c(176,183,173,191) # estatura de 4 homens
estatf <- c(157,152,174,166) # estatura de 4 mulheres
# Criando dois boxplot simples
boxplot(estatm, main="Boxplot de estatura de homens adultos (cm)",
xlab="Masculino", ylab="Estatura")boxplot(estatf, main="Boxplot de estatura de mulheres adultas (cm)",
xlab="Feminino", ylab="Estatura") Para juntar os dois gráficos, uma possibilidade é construir dois vetores de mesmo tamanho, um com os dados e outro que indique o sexo (Masculino ou Feminino):
estatm <- c(176,183,173,191)
estatf <- c(157,152,174,166)
estat <- c(estatm,estatf)
sexo <- rep(c("Masculino","Feminino"),each=4)
boxplot(estat~sexo, main="Boxplot de estatura de adultos (cm)", xlab="Sexo", ylab="Estatura (cm)")Experimente criar um arquivo .png:
png("estatura_homem_mulher.png")
boxplot(estat~sexo, main="Boxplot de estatura de adultos (cm)", xlab="Sexo", ylab="Estatura (cm)")
dev.off()A função png() tem vários outros parâmetros. Veja a documentação do R (com ?png): especialmente width e height são fundamentais. Há outros formatos gráficos disponíveis.
Você pode desviar a saída em tela para um arquivo texto. Experimente a função sink():
sink("estatura_homem_mulher.txt")
resumo <- summary(estat)
cat ("Estatura de todos os individuos:\n")
print (resumo)
sink()Precisa ser usada duas vezes: a primeira para abrir o arquivo e a segunda para fechá-lo - sink(), sem parâmetros - voltando a exibir as saídas na tela.
Na linha de métodos robustos é possível estimar o intervalo de confiança 95% por bootstrapping (https://en.wikipedia.org/wiki/Bootstrapping_(statistics)).
Por exemplo, para as estaturas dos homens deste exemplo:
rm <- replicate(1e4, mean(sample(estatm, replace=TRUE)))
plot(density(rm), main="Distribuicao das medias amostrais",
sub="(homens)", xlab="Media amostral", ylab="Densidade")qm <- quantile(rm, probs=c(0.025, 0.975))
cat (paste ("Intervalo de confianca 95% para os homens: [",qm[1],", ",qm[2], "]\n",sep=""))Intervalo de confianca 95% para os homens: [174.5, 187.25]
A função replicate(1e4, mean(sample(estatm, replace=TRUE))) é uma forma de fazer reamostragem (sample), de uma população hipotética com média (mean) igual à da amostra de estatura dos homens com reposição (replace=TRUE), repetindo (replicate) o processo 10.000 vezes (1e4).
O intervalo de confiança de 95% é obtido com a função quantile(), localizando os valores que deixam 2.5% da área sob as caudas esquerda e direita desta distribuição de probabilidades.
knitr::include_graphics("./image/William_Sealy_Gosset.jpg", dpi=80)Este teste foi inicialmente publicado por em 1908 por William Sealy Gosset sob o pseudônimo de Student, e posteriormente aprimorado por Ronald Fisher, que introduziu os graus de liberdade.
Em R, o teste \(t\) default é executado com a correção de Satterthwaite (erroneamente exibido como Welch), mas é possível forçar o teste clássico, adicionando o parâmetro var.equal:
Two Sample t-test
data: Sodium by Instructor
t = 0.76722, df = 38, p-value = 0.4477
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-67.59215 150.09215
sample estimates:
mean in group Brendon Small mean in group Coach McGuirk
1287.50 1246.25
O que muda são os graus de liberdade, aqui correspondendo a um número inteiro igual ao tamanho da amostra subtraído de duas unidades (uma para cada grupo).
Classicamente, havia várias condições para se poder executar o teste \(t\). No entanto, em suas versões robustas e com tamanho de amostra suficiente, podemos dispensar os testes prévios de homocedasticidade e normalidade. Quando estas condições não são atendidas, os pesquisadores podem optar pelos métodos equivalentes não paramétricos (e.g., Wilcoxon e Mann-Whitney).
Há respaldo na literatura especializada, sugerindo que atualmente esta pode não ser a melhor opção: